home *** CD-ROM | disk | FTP | other *** search
/ Atari Mega Archive 1 / Atari Mega Archive - Volume 1.iso / language / examples.zoo / numberth / lll.lsp < prev    next >
Lisp/Scheme  |  1991-10-22  |  7KB  |  137 lines
  1. ;; Lovász' Gitterbasis-Reduktions-Algorithmus
  2.  
  3. ;; Bruno Haible 15.8.1989
  4. ;; Nach
  5. ;;  A. K. Lenstra, H. W. Lenstra, L. Lovász: Factoring polynomials with
  6. ;;  rational coefficients. Math. Ann. 261 (1982), 515-534.
  7. ;; und
  8. ;;  Erich Kaltofen: Polynomial factorization 1982-1986, preprint.
  9.  
  10. (provide 'LLL)
  11.  
  12. ; Input: basis = Sequence von n linear unabhängigen Vektoren im Z^m
  13. ; Output: reduzierte Gitterbasis als Sequence von Vektoren im Z^m
  14. (defun LLL (basis)
  15.   (let ((n (length basis)))
  16.     (when (or (= n 0) (= n 1)) (return-from LLL basis))
  17.     (let ((m (length (elt basis 0))))
  18.       (labels ((scalprod (b1 b2) ; Skalarprodukt zweier Vektoren
  19.                  (let ((sum 0))
  20.                    (dotimes (i m) (incf sum (* (aref b1 i) (aref b2 i))))
  21.                    sum
  22.                ) )
  23.                (norm2 (b) ; Normquadrat eines Vektors
  24.                  (scalprod b b)
  25.               ))
  26.         ; Gram-Schmidt-Initialisierung:
  27.         (let ((b (map 'vector #'copy-seq basis))
  28.               (beta (make-array n))
  29.               (mu (make-array `(,n ,n))))
  30.           (do ((i 0 (1+ i))) ((= i n))
  31.             (do ((j 0 (1+ j))) ((= j i))
  32.               (setf (aref mu i j)
  33.                 (/ (- (scalprod (aref b i) (aref b j))
  34.                       (let ((sum 0))
  35.                         (do ((l 0 (1+ l))) ((= l j) sum)
  36.                           (incf sum (* (aref mu j l) (aref mu i l) (aref beta l)))
  37.                    )  ) )
  38.                    (aref beta j)
  39.             ) ) )
  40.             (setf (aref beta i)
  41.               (- (norm2 (aref b i))
  42.                  (let ((sum 0))
  43.                    (do ((l 0 (1+ l))) ((= l i) sum)
  44.                      (incf sum (* (expt (aref mu i l) 2) (aref beta l)))
  45.             ) )  ) )
  46.           )
  47.           ; Hier repräsentieren beta und mu die Koeffizienten der
  48.           ; Gram-Schmidt-Orthogonalisierung von b:
  49.           ; Für 0<=i<n ist b*[i] = b[i] - sum(j=0..i-1, mu[i,j] b*[j])
  50.           ;                mit  mu[i,j] = b[i] b*[j] / beta[j]
  51.           ;                und  beta[j] = b*[j] b*[j].
  52.           (do ((k 1)) (nil)
  53.             ; Hier ist das Teilstück b[0],...,b[k-1] der Basis reduziert,
  54.             ; d.h. für  0<=j<i<=k-1  ist  |mu[i,j]|<=1/2  und
  55.             ;      für  0<l<=k-1  ist  beta[l] >= beta[l-1]/2.
  56.             (if (or (= k 0) (>= (* 2 (aref beta k)) (aref beta (1- k))))
  57.               (progn
  58.                 ; Für  0<l<=k  ist  beta[l] >= beta[l-1]/2.
  59.                 (incf k) ; k erhöhen
  60.                 ; Für  0<l<k  ist  beta[l] >= beta[l-1]/2.
  61.                 (when (= k n) (return-from LLL b))
  62.                 ; Ersetze  b[k]  durch  b[k]-r*b[l]  für verschiedene
  63.                 ; r zu jedem l=0..k-1. Dabei bleiben alle b*[j] erhalten.
  64.                 (do ((l (1- k) (1- l))) ((< l 0))
  65.                   (let ((r (round (aref mu k l))))
  66.                     (unless (zerop r) ; bei r=0 verändert sich nichts
  67.                       ; (setf (aref b k) (- (aref b k) (* r (aref b l)))) ==
  68.                       (dotimes (j m)
  69.                         (decf (aref (aref b k) j) (* r (aref (aref b l) j)))
  70.                       )
  71.                       ; mu[k,.] adjustieren: Da jetzt gilt
  72.                       ; b*[k] = (b[k]+rb[l]) - sum(j=0..k-1, mu[k,j] b*[j])
  73.                       ;       = b[k] - sum(j=0..k-1, mu[k,j] b*[j])
  74.                       ;              + r b*[l] + r sum(j=0..l-1, mu[l,j] b*[j])
  75.                       ; muß mu[k,j] um r*mu[l,j] (für j<l) bzw. um r (für j=l)
  76.                       ; erniedrigt werden:
  77.                       (do ((j 0 (1+ j))) ((= j l))
  78.                         (decf (aref mu k j) (* r (aref mu l j)))
  79.                       )
  80.                       (decf (aref mu k l) r) ; dadurch wird |mu[k,l]|<=1/2
  81.               ) ) ) )
  82.               ; Vertausche b[k-1] und b[k]. Dabei sind auch beta[k-1] und
  83.               ; beta[k] zu vertauschen und mu[k-1,.], mu[k,.], mu[.,k-1],
  84.               ; mu[.,k] zu adjustieren.
  85.               (let* ((k-1 (1- k))
  86.                      (mu_k_k-1 (aref mu k k-1))
  87.                      (gamma_k-1 (+ (aref beta k) (* (expt mu_k_k-1 2) (aref beta k-1))))
  88.                      (h (/ (aref beta k) gamma_k-1))
  89.                      (nu_k_k-1 (/ (* mu_k_k-1 (aref beta k-1)) gamma_k-1)))
  90.                 (do ((j 0 (1+ j))) ((= j k-1))
  91.                   (rotatef (aref mu k-1 j) (aref mu k j))
  92.                 )
  93.                 (do ((i (1+ k) (1+ i))) ((= i n))
  94.                   (psetf (aref mu i k-1) (+ (* (aref mu i k-1) nu_k_k-1) (* (aref mu i k) h))
  95.                          (aref mu i k) (- (aref mu i k-1) (* (aref mu i k) mu_k_k-1))
  96.                 ) )
  97.                 (setf (aref beta k) (* (aref beta k-1) h))
  98.                 (setf (aref beta k-1) gamma_k-1)
  99.                 (setf (aref mu k k-1) nu_k_k-1)
  100.                 (rotatef (aref b k-1) (aref b k))
  101.                 (setq k k-1) ; Nun muß k verkleinert werden
  102.               )
  103.           ) )
  104. ) ) ) ) )
  105.  
  106. #|
  107. ; Bestimmt das Minimalpolynom (über Z) einer komplexen algebraischen Zahl x.
  108. ; Sie ist gegeben durch eine Funktion, die - mit einer ganzen Zahl k>=0 als
  109. ; Argument aufgerufen - eine komplexe Zahl liefert, die sich von 2^k*x sowohl
  110. ; im Realteil als auch im Imaginärteil um höchstens 1/2 unterscheidet.
  111. ; Dazu muß noch eine Gradabschätzung nach oben m und die L2-Norm |f| eines
  112. ; x annullierenden Polynoms f/=0 angegeben werden.
  113. ; Ergebnis ist das Minimalpolynom von x in der Form #(g0 ... gk)
  114. ; für g(x) = g0*x^0+...+gk*x^k == 0.
  115. (defun minpoly (fun m |f|)
  116.   (let* ((B0 (float |f| 0d0))
  117.          (B1 (* (float (/ (! (* 2 m)) (expt (! m) 2)) 0d0) B0))
  118.          (B2 (* (sqrt (+ m 2.0d0)) B1))
  119.          (B3 (* (sqrt (float (expt 2 m) 0d0)) B2))
  120.          (B4 (expt (* B1 B3) m))
  121.          (x0 (funcall fun 0)) ; grobe Näherung für x
  122.          (Betragsschranke ; Abschätzung für |x| nach oben, >=1
  123.              (max (sqrt (+ (expt (+ (abs (realpart x0)) 0.5d0) 2)
  124.                            (expt (+ (abs (imagpart x0)) 0.5d0) 2)
  125.                   )     )
  126.                   1.0d0
  127.          )   )
  128.          (B5 (* m (expt Betragsschranke (1- m)) B4))
  129.          (B6 (* (sqrt (+ m 3.0d0)) B3 B5))
  130.          (k (nth-value 1 (decode-float B6)) ; 2^k>B6
  131.          ; Brauche 2^k*x,...,2^k*x^m auf Genauigkeit 1/2.
  132.          (MM (* m (expt (+ Betragsschranke 0.25d0) (1- m)) (expt 2 (+ k 2))))
  133.          ; Brauche x auf Genauigkeit 1/MM.
  134.          (x (funcall fun (nth-value 1 (decode-float MM))))
  135. |#
  136.  
  137.